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ABSTRACT 

Wc develop an analytie model to calculate the rate at which galaxy disks are heated 
by dark matter substructures orbiting in their halos. The model takes into account 
the internal structure, mass function and accretion rate of satellites expected in the 
ACDM cosmology, as well as the growth of the disk by accretion and mergers, but 
it ignores resonant heating of the disk and the dynamical effects of spiral arms and 
bars. We calibrate this model against N-body simulations and demonstrate that it is 
able to reproduce the N-body heating rates to within a factor of 3 in the majority of 
cases. Our model gives the distribution of disk scale-heights for galaxies of different 
luminosities. For spiral galaxies, it predicts a median disk thickness of only 5% 
of the radial scale-length if substructure is the only source of heating. The median 
disk thickness increases to nearly 20% of the radial scale-length when heating due 
to gravitational scattering of stars by molecular clouds is also included. The latter 
value is close to the thickness estimated observationally for the disk of the Milky Way 
galaxy. The distribution of disk thickness predicted by the model is also consistent 
with a recent observational determination for sub-L* galaxies by Bizyaev & Mitronova 
(2002). Thus, the observed thickness of the stellar disks of spiral galaxies seems to be 
entirely compatible with the abundance of substructure in dark matter halos predicted 
by the standard A-dominated cold dark matter model of structure formation. In an 
fio = 1 universe, our best model of galaxy formation produces similar scale-heights, 
a consequence of the fact that similar amounts of substructure are accreted by halos 
during the lifetime of the disk in fio = 1 and JIq = 0.3, Aq = 0.7 cold dark matter 
cosmologies. 



1 INTRODUCTION 



A generic prediction of hierarchical models of structure for- 
mation, such as the cold dark matter (CDM) model, is that 
the dark matter halos of galaxies and clusters should con- 
tain large amounts of substructure, in the form of small, 
gravitationally bound subhalos orbiting within the larger 
potential. This substructure arises because large halos are 
built up by mergers of smaller halos whose tidally-stripped 
remnants can survive in favourable conditions. Recently, it 
has been claimed that the CDM model predicts an order of 
magnitude too many subhalos around the Milky Way galaxy, 
compared to what is inferred from the number of satellite 
galajcies (Klypin et al. 1999; Moore etal 1999). Several au- 
thors have now pointed out that this apparent discrepancy is 
readily explained if some process (such as the heating of the 
intergalactic medium (IGM) during reionization) is efficient 
at suppressing the formation of galaxies in most of these 
subhalos (Bullock, Kravtsov & Weinberg 2000; Benson et 
al. 2002b; Somerville 2002). In this picture, galaxy halos 
should be filled with many small subhalos containing neg- 



figible amounts of luminous material. A strong test of this 
idea is possible by searching for gravitational signatures of 
subhalos, thus bypassing the problem of relating subhalos 
to the visible material in satellite galaxies. 

The most direct probe of substructure in dark matter 
halos is gravitational microlensing. Its properties are reason- 
ably well understood theoretically (Mao & Schneider 1998; 
Metcalf & Madau 2001; Chiba 2002; Dalai & Kochanek 
2002a; Dalai & Kochanek 2002b). Although the interpre- 
tation of the current datasets remains controversial in some 
cases, the observed microlensing rates appear to be consis- 
tent with the abundance of substructure predicted by CDM. 

An alternative constraint on the amount of substruc- 
ture in halos may be obtained by considering the thickness 
of the stellar disks of galaxies. Subhalos on orbits that pass 
through or near to a galactic disk perturb it gravitationally 
and deposit energy into it, gradually heating the disk and 
increasing its scale-height. Since there are other mechanisms 
that also heat stellar disks (but with uncertain efficiency), 
the observed thickness of galactic disks sets an upper limit 
on the abundance of such substructure. The heating of galac- 
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tic disks by infalling satellites was invoked as a constraint 
on models of structure formation by Toth & Ostriker (1992; 
hereafter TO). They calculated this effect semi-analytically, 
and concluded that the thinness of the Milky Way's disk 
is inconsistent with the hierarchical build-up of galaxies in 
a high density (f2o = 1) CDM universe. This conclusion 
was disputed by Navarro, Frenk & White (1994) whose 
cosmological simulations showed that many of the satellites 
that are incorporated into a galactic dark halo do not ac- 
tually merge with the central galaxy. Subsequent numerical 
simulations of mergers of single satellites with larger disk 
galaxies (e.g. Huang & Carlberg 1997; Velazquez & White 
1999) indicated that TO's analytical estimates of the heat- 
ing rate were somewhat too high, weakening their constraint 
on structure formation models. More recently. Font et al. 
(2001) have numerically simulated the heating of disks by 
the ensemble of subhalos predicted to exist within dark ha- 
los in the CDM model. Their simulations of Milky Way-like 
galaxies only set an upper limit to the rate of disk heating 
by satellites, because of numerical effects, but they conclude 
that this is less than the total disk heating rate that is in- 
ferred observationally for the Solar neighbourhood. They ar- 
gue that the heating rates are low because the most massive 
satellites, which are the ones that cause the most heating, 
are few in number and because few satellites penetrate the 
inner regions of the disk. Although their conclusions agree 
with those of Navarro, Frenk & White (1994), they are lim- 
ited by the fact that they simulated only two realizations of 
the halo substructure. 

In this paper, we develop a new semi-analytical model 
of disk heating by halo substructure. Our calculation builds 
upon earlier semi-analytical modeling of galaxy formation 
within the framework of CDM cosmology (Cole et al. 2000), 
and on recently developed analytical models of the evolution 
of satellites within dark matter halos (Taylor & Babul 2001; 
Benson et al. 2002a; Taffoni et al. 2003; Taylor & Babul 
2003). The rate at which satellite halos of different masses 
are incorporated into the main halo is given by the galaxy 
formation model. The satellite model then predicts how the 
masses, radii and orbits of subhalos evolve due to dynamical 
friction and tidal stripping by the halo, disk and bulge of the 
host galaxy. In this paper, we add a calculation of how much 
of the orbital energy of the satellites that is lost by dynam- 
ical friction goes into increasing the thickness and vertical 
motions of the galactic disk. The interaction between the 
satellite and the disk is modeled in a simplified way, ignor- 
ing details such as resonant interactions and the possible 
role of spiral features and bars. We test and calibrate our 
analytical model of satellite evolution against a new set of 
high-resolution N-body simulations of single satellites merg- 
ing with disk galaxies. We find (as has also been shown by 
Taylor & Babul 2001 and Taffoni et al. 2003) that such an 
analytical model is able to reproduce well the behaviour seen 
in the N-body simulations. We measure the disk heating in 
the same simulations, and find that it is quite well repro- 
duced by our analytical model. We then apply this model 
of heating by satellites within the framework of our semi- 
analytical model of galaxy formation, in order to predict 



the distribution of scale-heights for disk galaxies of different 
luminosities. 

Both the N-body and semi- analytical approaches have 
advantages and disadvantages when applied to this problem. 
N-body simulations fully account for the non-linear inter- 
action of substructure and disk (e.g. for the excitation of 
global modes such as warps and bars in the disk). However, 
they are limited by resolution and artificial numerical heat- 
ing and, because of computational cost, they are limited to 
few (two, in the case of the best cosmological simulations of 
disk heating to date, by Font et al. ). The semi-analytical 
approach has the advantage that it is not limited by resolu- 
tion or artificial heating, and it allows the calculation of a 
large number of realizations. Since heating by substructure 
is a highly stochastic process, it is important to account for 
the galaxy-to-galaxy variation in the heating rate by calcu- 
lating a large number of realizations. At present, this is only 
possible with the semi-analytical approach. 

The remainder of this paper is arranged as follows. In 
§2 we describe our analytical model for disk heating by sub- 
halos and for the evolution of the disk scale-height. In §3 we 
calibrate and test our analytical model against numerical 
simulations of single satellite-disk mergers. In §4 we present 
our predictions for the distribution of scale-heights of disk 
galaxies in the CDM model, and compare with observational 
data for the Milky Way and for other galaxies. Finally, in §5 
we present our conclusions. Appendices detail derivations 
of various formulae related to dynamical friction and disk 
energies and present convergence tests for the N-body sim- 
ulations. 



2 MODEL 

2.1 Evolution of Satellites and their Orbits 

We calculate the evolution of the masses, radii and orbits 
of satellites using a development of the model presented in 
Benson et al. (2002a; hereafter Paper I). That work, in turn, 
was based on the satellite evolution model of Taylor & Babul 
(2001). We summarize here the main features of our model. 
The growth of the main halo is described by a merger his- 
tory tree which is calculated by a Monte Carlo method (Cole 
et al. 2000). When smaller halos (in general containing one 
or more visible galaxies) merge with the main halo, they 
become satellite halos. The satellite halos are given initial 
orbits which start close to the virial radius but have a range 
of eccentricities consistent with the distribution seen in the 
N-body simulations of Ghigna et al. (1998.). The satellite 
orbits are followed in the potential of the host system; they 
evolve due to dynamical friction against the dark halo, disk 
and bulge of the main galaxy. At the same time, the satel- 
lites lose mass by tidal stripping, both "static" tidal limita- 
tion and tidal shocking. As a satellite is tidally stripped, its 
radius and internal structure also change. 

We have made a few improvements to our satellite orbit 
model from that presented in Paper I. These are described 
in Appendix A. 
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2.2 Disk Heating 

2.2.1 Rate of Heating 

We now wish to calculate the rate at which a satellite halo 
heats the disk of the galaxy in its host halo. The satellite ex- 
periences dynamical friction against the disk, and the energy 
lost from the orbital motion of the satellite by this mecha- 
nism goes into increasing the energy of the disk. Working in 
the frame in which the centre of mass of the central galaxy 
and its halo are at rest, the satellite injects energy into the 
disk at a rate 

P = -Fdf.disk • Vsat, (1) 

where Fdf,disk is the dynamical friction force exerted by the 
disk and Vaat is the velocity of satellite. (Note that while we 
typically expect the satellite to lose energy to the disk, it 
is possible for the satellite to gain energy from the ordered 
motions of the disk if Fdf.disk • Vsat > 0. This occurs because 
the dynamical friction force depends, in our approximation, 
on the relative velocity vector of the satellite and the local 
disk stars. If the local disk velocity is sufficiently large, the 
relative velocity vector may point in the opposite direction 
to the satellite velocity vector, resulting in a transfer of en- 
ergy from the disk to the satellite.) This energy is initially 
injected in the form of kinetic energy, but it is subsequently 
mixed between kinetic and potential energies by the motions 
of the stars. We are interested in the increase in the vertical 
energy of the disk, which is given by 

Ez = — tzFdf.disk • Vsat, (2) 

We derive an expression for the efficiency factor < 1 in 
Appendix B3 by considering the increases in the vertical and 
horizontal components of the velocity dispersion of the stars 
responsible for dynamical friction during scattering events. 
This expression (B29) depends only on the Coulomb log- 
arithm. In A, and the angle 6vg between the disk-satellite 
relative velocity and the z-axis. We then simply integrate 
Ez along the satellite orbit to determine the net increase in 
the disk's vertical energy. 

Figure 1 shows how depends on the angle dy^ for 
a few representative values of A. Note that ez = | when 
cos dvg = , independently of A. For small A the effi- 
ciency is greatest when Oyg = 0° (approaching unity as A 
approaches zero) and smallest for ^Vq = 90° (approaching 
zero as A approaches zero). For large A the trend is reversed, 
with ez being smallest at Oy^ = 0° (approaching zero as A 
approaches infinity) and largest for Oy^ = 90° (approaching 
i as A approaches infinity). The transition between these 
two regimes occurs for A ~ 3.975, for which ez is indepen- 
dent of 9yg . 

We can understand the behaviour of ez in simple terms. 
For example, for 9yg — 0, the efficiency drops to zero as 
A becomes large. In this case, vertical motions in the disk 
are parallel to the relative velocity vector of the satellite 
and the disk stars. Consequently, only the AT^-n|| term (see 
eqn. B15) contributes to increasing the energy in these ver- 
tical motions. As A, and hence the maximum impact param- 
eter, b, increases, energy transfer from the satellite becomes 
dominated by large b scatterings. For large impact parame- 
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Figure 1. The cflicioncy of energy transfer to vertical motions in 
the disk as a function of the angle between the disk-satellite rel- 
ative velocity and the z-axis, dy^ . Results are plotted for several 
values of A as indicated in the figure. 



ters, the increase in velocities (and hence energies) perpen- 
dicular to the satellite motion dominates over that parallel 
to the motion, since Al/mii oc b~^ while AKn± k b~^ (see 
eqn. B16). Consequently, the efficiency of transfer to vertical 
motions in the disk drops to zero as A becomes large. 

The reversal of the trend of ez with 6'vo at A ~ 3.975 
is also simple to understand. For larger A, energy transfer 
is predominantly into motions perpendicular to the motion 
of the satellite (as discussed above). Thus, the efficiency 
of energy transfer to motions in the vertical direction is 
greatest when the satellite moves perpendicular to that di- 
rection {Oyg = 90°). For smaller A, energy transfer occurs 
mostly into the parallel direction, and so ez is maximized for 
^Vo = 0°. For A ~ 3.975 energy transfers into perpendicular 
and parallel directions are equal and so ez is constant. 

When cos 6yg = , the energy transferred to vertical 
motions is always one third of the increase in energy parallel 
to the satellite velocity, plus one third of the increase in the 
energies in the two directions perpendicular to the satellite 
velocity. Thus, this energy is always exactly one third of the 
total energy transferred from the satellite and hence ez = | 
independently of A. 

It is worth considering at this point some of the sim- 
plifications which go into our dynamical model of disk heat- 
ing. Dynamical friction is treated using Chandrasekhar's ap- 
proximation which is clearly not strictly applicable to our 
halo-plus-disk system. While this approximation has been 
shown to be a reasonable one for dark matter halos (Wein- 
berg 1986; Bontekoe & van Albada 1987; Core, Muzzio 
& Vergne 1997; van den Bosch 1999; Velazquez & White 
1999; Colpi, Mayer & Governato 1999), its validity when 
disks are included is less clear. Importantly, this approxima- 
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Table 1. Properties of the galaxy and satellite models used in the N-body simulations. The first column specifies the component in 
question. The second column gives the density profile, cither in spherical coordinates (r), or cylindrical polar coordinates {R,z). The 
values of the parameters of each profile are listed in column 3. The satellites are all described by King models (King f966). For these, 
we specify the core radius, t^, and the concentration c = \og^Qrt/r^ where rt is the tidal radius of the satellite. The final column lists 
the number of particles used to represent each component in the standard case. 



Component Density Profile 



Parameters 



Number of Particles 



Halo 



Disk 



Bulge 



Ph(r-) 



M, 



Ph{r) 



My. 



27V r(a + r)'^ 

Satellite SI King Model 
Satellite S2 King Model 
Satellite S3 King Model 



Mh = 7.84 X 10" M0 
7 = 3.5kpc 
r-cut = 84kpc 
a = 1.076 

^ exp(-if,/i?.d)soch2(2/Hd) Md = 5.6 X 10^°Mq 

i?d = 3.5kpc 
Hd = 700pc 
Mb = 1.87 X IO^OMq 
a = 525pc 

Ms = 5.60 X ID^Mo 
Tc = Ikpc 
c = 0.8 

Ma = 5.60 X lO^M© 
rc = 500pc 
c = 1.1 

Ma = 1.12 X IO^OMq 
rc = 875pc 
c = 1.0 



687008 

163840 

16384 
32768 

32768 

32768 



Table 2. Properties and initial orbital parameters of the satel- 
lites in the N-body simulations. Column 2 specifies the satellite 
model used (as defined in Table 1). Column 3 lists 9i, the angle 
between the initial angular momentum vector of the satellite and 
that of the disk. Column 4 lists the circularity of the satellite's 
initial orbit, ej, while column 5 lists the initial radial position of 
the satellite (which is the apocentre of its orbit), ra. Column 6 
specifies whether the simulation contains a disk or not (note that 
6i is undefined for diskless simulations G2Sxx). 



Model Satellite 



ra/kpc Disk? 



GlSl 


SI 


45° 


0.33 


59.0 




G1S2 


SI 


0° 


0.55 


55.0 


V 


G1S3 


SI 


45° 


0.55 


55.0 


V 


G1S4 


SI 


90° 


0.55 


55.0 


V 


G1S5 


SI 


135° 


0.55 


55.0 


V 


G1S6 


SI 


180° 


0.55 


55.0 


V 


G1S7 


SI 


0° 


0.82 


46.5 


V 


G1S8 


SI 


45° 


0.82 


46.5 


V 


G1S9 


S2 


0° 


0.55 


55.0 


V 


GISIO 


S2 


45° 


0.55 


55.0 


V 


GlSll 


S2 


90° 


0.55 


55.0 


V 


G1S12 


S2 


135° 


0.55 


55.0 


V 


G1S13 


S2 


180° 


0.55 


55.0 


V 


G1S14 


S3 


45° 


0.55 


55.0 


V 


G1S15 


S3 


135° 


0.55 


55.0 


V 


G2S1 


SI 


N/A 


0.33 


59.0 


X 


G2S2 


SI 


N/A 


0.55 


55.0 


X 


G2S7 


SI 


N/A 


0.82 


46.5 


X 


G2S9 


S2 


N/A 


0.55 


55.0 


X 


G2S14 


S3 


N/A 


0.55 


55.0 


X 



tion ignores any possible resonant interaction between the 
satellite and the disk. 



A further simplification of our model is that the phase 
space distributions of halo dark matter and disk stars are 
assumed to be fixed, with the exception that the disk vertical 
velocity dispersion and density profile are allowed to change 
with time. (We further assume that the vertical motions 
of stars in the disk do not couple to radial and azimuthal 
motions, which will be approximately true provided that 
the disk remains thin.) In reality, all three components of 
the disk velocity dispersion will be affected by substructure 
heating. However, the changes in the radial and azimuthal 
velocity dispersions have only a small effect on the overall 
structure of the disk in the majority of cases. Thus, our 
approach should be a reasonable first approximation. 

A final, important simplification is that we ignore some 
possible interactions between the disk and the dark matter 
halo, e.g. those driven by non-axisymmetric structures such 
as bars or warps in the disk. This complex set of interactions 
could, in principle, result in energy initially transferred from 
the satellite to the disk finding its way into the halo dark 
matter. The efficiency with which this happens will clearly 
depend upon the frequency with which substructure excites 
bars and other global modes in the disk and is therefore 
beyond the scope of our current calculations. 

Given these simplifications it is important to test our 
analytic calculations against N-body simulations of the disk 
heating process. We perform such tests in §3. 
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Table 3. Comparison of results from the analytic and N-body 
calculations of satellite evolution. Column 1 lists the model num- 
ber; column 2 lists the type of calculation (analytic or N-body); 
column 3 gives the change in the disk vertical kinetic energy gen- 
erated by the satellite at the end of the simulation, both in abso- 
lute units and as a percentage of the initial disk vertical kinetic 
energy (values in parentheses). Where the analytic and N-body 
estimates of the disk energy differ by more than a factor of two, 
we show the values in bold type. Columns 4 and 5 list the times 
at which the satellite reaches 50% and 10% of its initial mass 
respectively. 



Model 


Type 


10 


AT^(4Gyr) 


tso/Gyr 


tio/Gyr 


GlSl 


Analytic 





133 


(25%) 


1.24 


2.01 


GlSl 


N-body 





101 


(18%) 


1.38 


2.16 


G1S2 


Aiiajlytic 





257 


(49%) 


1.51 


2.89 


G1S2 


N-body 





342 


(60%) 


1.69 


2.08 


G1S3 


Anajlytric 





184 


(35%) 


2.20 


3.19 


G1S3 


N-body 





238 


(42%) 


2.25 


3.13 


G1S4 


Anailytric 





029 


(5%) 


2.32 


3.73 


G1S4 


N-body 





101 


(18%) 


2.54 


3.53 


G1S5 


Anajlytric 





032 


(6%) 


2.24 


3.77 


G1S5 


N-body 





057 


(10%) 


2.46 


3.55 


G1S6 


Ann l^/tir" 





177 


(33%) 


1.94 


3.00 


G1S6 


N-body 





090 


(16%) 


2.21 


3.14 


G1S7 


Analytic 





056 


(11%) 


2.70 


> 4.00 


G1S7 


N-body 





443 


(77%) 


2.22 


2.28 


G1S8 


Analytic 





096 


(18%) 


3.23 


> 4.00 


G1S8 


N-body 





324 


(57%) 


3.72 


4.03 


G1S9 


Analytic 





272 


(51%) 


1.92 


3.30 


G1S9 


N-body 





307 


(54%) 


1.82 


1.88 


GISIO 


Analytic 





244 


(46%) 


2.53 


3.28 


GISIO 


N-body 





588 


(103%) 


3.01 


3.35 


GlSll 


Analytic 





114 


(22%) 


2.86 


3.59 


GlSll 


N-body 





363 


(63%) 


3.16 


3.80 


G1S12 


Analytic 





131 


(25%) 


2.89 


3.58 


G1S12 


N-body 





229 


(40%) 


3.26 


4.06 


G1S13 


Analytic 





507 


(96%) 


2.32 


2.62 


G1S13 


N-body 





350 


(61%) 


2.87 


3.28 


G1S14 


Analytic 





521 


(98%) 


1.58 


2.52 


G1S14 


N-body 





873 


(153%) 


1.62 


1.90 


G1S15 
G1S15 


Analytic 
N-body 






438 
374 


(83%) 
(65%) 


1.78 
1.80 


2.09 
2.18 



2.2.2 Disk Scale-Height and Vertical Energy 

Having calculated the energy deposited into vertical motions 
of disk stars, we now wish to calculate the resulting scale- 
height of the disk. We work throughout in the thin disk ap- 
proximation, in which the vertical extent of the disk is always 
assumed to be small compared to its radial extent, and the 



non-circular velocities are assumed to be small compared to 
the circular velocity. In this approximation, the disk can be 
treated as being locally plane-parallel, with the consequence 
that the vertical motions separate from the motions in the 
plane, and there is a well-defined vertical energy which (in 
the absence of perturbations by satellites or other objects) 
is conserved both for individual stars and for the disk as a 
whole. The vertical energy given to a star by an encounter 
with a satellite is initially in the form of vertical kinetic en- 
ergy, but the orbital motion of the star subsequently mixes 
this between vertical kinetic and potential energies, while 
keeping the sum of the kinetic and potential energies con- 
stant. In the thin disk approximation, the total vertical en- 
ergy per unit area of the disk, e^, can be written as (TO) 



fiz = iz + Wdd + Wdh, 



(3) 



where all quantities are surface energy densities, tz is the 
disk vertical kinetic energy, Wdd is the disk self-gravitational 
energy and Wdh is the gravitational energy due to the 
disk/halo interaction. The vertical energy Cz is defined to 
be zero in a state where the disk has zero thickness and 
zero vertical velocities. Expressions for tz, Wdd and Wdh are 
derived in Appendix C. Following TO, we assume virial equi- 
librium and find 



2tz — Wdd - 

and so 
3 



2wdh = 0, 



-Wdd + 2wdh. 



(4) 



(5) 



The density of our model disks in the vertical direction is 
proportional to sech^z/Hd. For this density profile we find 
from eqn.(5) (TO) 

Q „2 



-■KGEi(R)Hd 



-E,(i?)i/,— 



(6) 



where R is radius in the disk plane, Ed(i?) is the disk surface 
mass density, and M\^{FVj is the mass in the (spherical) halo 
plus bulge within radius R. Since the vertical kinetic energy 
per unit area is tz = ISdtr^, we also find from eqn.(4) 

2 

ol = -nGT.d{R)Hd + —GMi,{R)Hl/R' (7) 

This expression is used to calculate the vertical velocity dis- 
persion at each radius from the scale-height Hd. * . 

To relate the radially-dependent vertical energy per unit 
area to the global total vertical energy, we make the as- 
sumption that the disk scale-height is constant with radius, 
since this is observed to be a good approximation for real 
galaxies (e.g. de Grijs & Peletier 1997). We can then inte- 
grate eqn.(6) over the whole disk to find the total vertical 
energy. Using Ed ~ (Md/iTrRd) exp{—R/ Rd) for an expo- 
nential disk of radial scale-length Rd we find 







Jo 


.Vd. 



exp(- 



-da;, (8) 



—MdVd'h+—MdVd'h' 



Note that here we differ slightly from VW by including the 
contribution of the halo gravity to the disk vertical velocity dis- 
persion. This is typically a small, although not negligible, contri- 
bution over the bulk of the disk. 
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where the fractional scale-height h = H^/Rd, = 
GAfd/-Rd and = GMb(R)/R. Integrating eqn.(7) gives 
a similar expression for the total vertical kinetic energy Tz. 
Once the total vertical energy is known, the above equa- 
tion is easily solved for h and hence _ffd- 



2.2.3 Local vs. Global Heating 

In §2.2.2 we made the assumption that the energy deposited 
in the disk by satellites was distributed throughout the disk 
in such a way as to produce a scale-height that was indepen- 
dent of radius. However, the increase in energy per unit mass 
caused by a satellite passing through or near the disk will 
initially be greatest close to the point of impact. Since satel- 
lite encounters frequently trigger global modes of the disk 
it is not implausible that this energy quickly becomes redis- 
tributed throughout the disk. However, it is interesting to 
consider the opposite extreme in which energy is deposited 
at the position of the satellite and remains there. We refer 
to these two extremes as "global" and "local" heating. To 
study local heating we accumulate the energy deposited by 
satellites in a narrow annulus of the disk (in practice we use 
a Gaussian window function), centred on the disk half-mass 
radius. We then assume that the specific energy of disk ma- 
terial is proportional to the same window function and use 
the relations of §2.2.2 to compute the resulting scale-height 
at the half-mass radius. 

Observations of real galactic disks (de Grijs & Peletier 
1997) indicate that the scale-height is reasonably constant 
with radius, at least for late-type galaxies. For this reason 
we prefer the global heating assumption, but consider local 
heating also as an interesting comparison. 



2.2.4 Further Aspects 

Below we detail how we deal with energy^ deposited in a 
gaseous disk and how we treat galaxy mergers, gas accretion 
and star formation. 

Gas in Galaxy Disks: Disks in our model in general con- 
sist of both stars and gas. The gas is assumed to be in an 
infinitely thin layer with zero velocity dispersion in the disk 
midplane. We include the contribution of the gas to the disk 
gravitational potential and when computing the disk scale- 
height. With our choice of zero-points for the energy, the ver- 
tical kinetic energy of the gas and also its self-gravitational 
energy are both zero (because it is at 2; = 0), but the gas 
contributes to the total energy per unit area of the disk 
ez through the gravitational interaction energy between the 
gas and stars (see TO for more details) . We assume that gas 
and stars in the disk are heated at the same rate per unit 
mass, but that the gas dissipates this energy rapidly, so that 
energy deposited in the gas is effectively lost. 

Adiabatic Heating due to Gas Accretion: Gas accreted 



t For convenience, we use the expression "energy" to imply "disk 
vertical energy" from here on, unless explicitly stated otherwise. 



onto the disk is assumed to initially have zero energy. How- 
ever, the growth of the disk surface density causes gravita- 
tional compression in the vertical direction, which tends to 
increase the vertical energies of disk stars. We follow TO 
and assume that gas infall occurs adiabatically, adopting 
their equation (3.12) to describe the change in energy of the 
disk stars due to adiabatic heating. In our model, gas can 
also be lost from the disk due to feedback processes, result- 
ing in a decrease in the energies of stars. We account for this 
process in the same way as for the adiabatic heating, simply 
changing the sign of the effect. We find that these are minor 
effects that have little impact on the predicted scale-heights 
of disk galaxies. 

Star Formation: When gas turns into stars, we assume 
that these stars start out with zero energy, but then rapidly 
mix with the pre-existing stellar population, conserving the 
total disk vertical energy. 

Galaxy Mergers: In a major merger all disks are de- 
stroyed, and so we zero the energy of the resulting galaxy. 
In minor mergers, stars from the satellite galaxy disk and 
bulge are added to the bulge of the central galaxy. In the 
merger, the energy of the satellite disk is lost, while that 
of the central galaxy disk is unchanged, unless the infalling 
satellite contains gaseous material, which will contribute to 
the adiabatic heating of the central galaxy disk. 



2.2.5 Heating of Disks by Scattering by Clouds 

Substructure in the halo is not the only source of heating 
for disks. Two other plausible mechanisms for disk heating 
are gravitational scattering of stars by massive gas clouds 
(Spitzer & Schwarzchild 1953; Lacey 1984) and scattering 
of stars by spiral arms (Carlberg & Sellwood 1985). The 
latter mechanism is inefficient at producing any heating in 
the vertical direction, so we will focus on the first mecha- 
nism. Lacey (1984) derived analytical expressions for the 
rate at which scattering by clouds increases the vertical and 
horizontal epicyclic energies of the stars. In general, these 
expressions depend on the radial and vertical disk velocity 
dispersions, ctr and a^, but, acting by themselves, the clouds 
tend to drive the ratio az/an to an equilibrium value. We 
calculate the rate of increase of vertical energy per unit mass 
for the stars, Ez, using Lacey's eqn.(39), evaluated for the 
equilibrium (Tz/o-r and in the limit in which the scale-height 
of the stars is larger than that of the clouds. This gives 



/dez 
V dt 



2 a^{l3)Ks(P) 



(9) 



where Ec is the surface density in clouds, Mc is the cloud 
mass. In Ac is the Coulomb logarithm for scattering of stars 
by clouds and v is the vertical epicyclic frequency. as(/?) and 
Ks{l3) axe functions of /3 = 2il/K which are tabulated by 
Lacey, Q being the angular velocity for circular orbits and n 
the radial epicyclic frequency. We obtain the total contribu- 
tion of scattering by clouds to the increase of vertical energy 
by integrating eqn.(9) over radius: 



z. clouds — 



4w) 

\ at / cloud 



2ttR dR 



(10) 
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Numerical simulations of heating by clouds agree fairly 
well with the velocity dependence predicted analytically, 
(d(T^/dt oc o-~^), but have given somewhat conflicting results 
about the amplitude of the effect; Villumsen (1985) found 
heating rates da^ /dt at a given a about 6 times lower than 
the analytical prediction, while Hanninen & Flynn (2002) 
found rates 3-8 times higher. 

Our galaxy formation model predicts the total mass of 
gas in the disk of each galaxy as a function of time. We as- 
sume that the gas is distributed radially in the same way 
as the stars, with a constant fraction being in the form of 
giant molecular clouds. For our standard case we will as- 
sume that 25% of the gaseous mass of the disk is in clouds 
(Granato et al. 2000), that they have mass- weighted mean 
mass of A4c = 6.6 x 10^ M0 (Lacey 1984) and typical radius 
Oc = 16pc (Granato et al. 2000), and that /? — 1.5. For 
each model galaxy, we integrate the heating due to scatter- 
ing from molecular clouds over each timestep in the calcula- 
tions, and add this energy change to that which arises from 
interactions with satellites. 



3 CALIBRATION USING N-BODY 
SIMULATIONS 

As has been noted by several authors, the amount of heat- 
ing caused by a satellite is difficult to determine analytically 
since some of the energy may drive global perturbations (e.g. 
warps) in the disk, and satellites may trigger bar instabili- 
ties leading to an enhanced heating rate. Furthermore, our 
approach to dynamical friction in the disk follows the meth- 
ods of Chandrasekhar (e.g. Binney & Tremaine 1987, sec- 
tion 7.1), which assume that each particle interacts with the 
satellite only once. If the satellite orbital period is close to 
the rotation period of the disk (or to some other resonance 
of the disk orbits), this assumption fails. Instead, a single 
particle may interact with the satellite multiple times on 
consecutive orbits. This problem should therefore ideally be 
approached in terms of resonant interactions between satel- 
lite and disk (Goldreich & Tremaine 1979; Donner & Sun- 
delius 1993; Wahde, Donner & Sundelius 1996; Weinberg 
& Katz 2002). We retain the Chandrasekhar methods for 
their simplicity, and show that they provide a reasonable 
approximation to the dynamical friction due to disks in the 
regimes of interest. 



3.1 N-body Simulations 

We begin by testing and calibrating our analytic calcula- 
tions against numerical simulations of disk heating. In prin- 
ciple, the simulations of VW are ideal for this purpose. How- 
ever, the central densities and velocity dispersions of the 
King model satellites given by VW are too low to be con- 
sistent with their assumed concentration parameters. Thus, 
the satellites seem to more weakly bound than the authors 
intended. It is unclear a priori how this would affect the 
results and we have therefore decided to repeat their calcu- 
lations. This has two other advantages: 



• We can repeat each simulation without the disk compo- 
nent, allowing us to constrain separately the contributions 
of the halo and disk to the dynamical friction experienced 
by the satellite. 

• We can perform convergence tests by increasing the 
number of particles in the simulation in order to ensure that 
disk heating is being estimated accurately. 

We carry out the same set of simulations as VW. Brieffy, 
each simulation consists of a galaxy containing a bulge, 
disk and dark matter halo, plus a satellite object. Density 
profiles and the number of particles used for each compo- 
nent are listed in Table 1, while other details of each sim- 
ulation (type of satellite used, initial satellite orbital pa- 
rameters and whether or not a disk is included) are listed 
in Table 2. Initial conditions are created using the tech- 
niques of Hernquist (1993). The galaxy and satellite are 
then evolved separately, as described by VW, using the 
GADGET code (Springel, Yoshida & White 2001) to allow 
them to reach equilibrium. We employ gadget's new cell- 
opening criterion for tree walks (TypeOf OpeningCriteria= 
1) with an accuracy of ErrTolForceAcc= 0.001, together 
with TypeOf TimestepCriterion= 1 with ErrTolVelScale= 
10.0. GADGET uses adaptive time-stepping. We impose 
no minimum timestep size, but impose a maximum size 
of MaxSizeTimestep= 0.01 (in GADGET'S default internal 
units). All particles in the simulation are given a softening 
length of O.llOkpc. With these choices, energy is conserved 
to better than 1% throughout the simulations. The two sets 
of initial conditions are then superimposed and evolved for 
4Gyr. 

The simulations are labelled GlSl to G1S15 as in VW. 
We also perform a simulation with no satellite, GISO, to 
measure the two-body heating rate in the disk. We repeat 
each simulation without a disk component, labelling these 
G2S1 to C2S15 (note that in the absence of the disk, only 
models G2S1, G2S2, G2S7, G2S9 and G2S14 are different). 
We also repeated all of these calculations with one half and 
one quarter the number of particles, in order to test how 
well the results have converged. The convergence tests are 
described in Appendix D. They indicate that the conver- 
gence is good for the evolution of the mass and orbit of the 
satellite, and adequate for the increase in the vertical energy 
of the disk. Unless otherwise noted, we show results from the 
highest resolution simulations. 

Each simulation output is analyzed in order to deter- 
mine the position, velocity and mass of the satellite (com- 
puted for those particles which remain bound to the satel- 
lite), and the vertical kinetic energy of the disk. We deter- 
mine which particles are bound to the satellite using the 
following algorithm: 

(i) Begin by considering all the satellite particles that 
were bound to the satellite at the previous timestep (or sim- 
ply all satellite particles for the first timestep). 

(ii) Compute the mean position and velocity, and the 
mass of the satellite from these particles. 

(iii) For each particle in this set, determine if it is gravi- 
tationally bound to the other particles in the set. 

(iv) Retain only those particles which are bound and go 
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back to step (ii). Repeat until the mass of the sateUite has 
converged. 

To determine the vertical kinetic energy of the disk, 
Tz, at each output time, we first locate the centre of mass 
of the disk and its mean velocity. (Since the satellite mass 
is comparable to that of the disk, the disk moves around 
significantly as the satellite passes by.) We then rotate the 
system to the frame defined by the principal axes of the disk 
inertia tensor, and sum the kinetic energies of particles in 
the direction defined by the shortest axis (which corresponds 
to the z-axis for an untilted disk). This rotation is necessary 
because the disk can become tilted through its interaction 
with the satellite (as also noted by VW). In the original 
frame (i.e. without rotation), purely circular motions in a 
tilted disk appear as vertical energy. 

One final step is necessary in order to obtain the in- 
crease in the disk vertical energy due to the interaction 
with the satellite. Even in the absence of a satellite, the 
vertical kinetic energy of the disk increases as the simu- 
lation proceeds due to numerical relaxation (mainly two- 
body scattering), from 0.57 x lO^Mokm^s^^ at f = to 
0.64 X lO^''M0km^s~^ at t = 4Gyr, for our standard particle 
number. This increase of 0.07 x lO^^MQkm'^s"'^ due to two- 
body relaxation is comparable to the heating by the satellite 
in many of the cases considered. Therefore, to obtain the in- 
crease in vertical energy due to the satellite at time t, which 
we denote as ATz it) , we subtract off the energy of the un- 
perturbed disk (from model GISO) at the same time t. Based 
on runs of model GISO with different random number seeds 
but the same number of particles, the increase in Tz due to 
numerical relaxation is determined to an accuracy of better 
than 0.005 x lO^Mekm^s"^ in the standard case, so the 
uncertainty in ATz introduced by the subtraction is small 
compared to ATz itself. 

3.2 Comparison with Analytic Calculations 

To test our analytical model of satellite orbital evolution and 
constrain its parameters, we adapt the analytical model so as 
to mimic the set-up of each N-body simulation. Thus, we as- 
sume density profiles for host and satellite systems identical 
to those of the N-body simulations. Taylor & Babul (2001) 
compared their model of satellite galaxy orbital evolution to 
the orbital radii and satellite masses as a function of time in 
the simulations of VW, finding generally good agreement. 
We repeat their analysis here, using our own model of satel- 
lite dynamics, extended to include the calculation of disk 
heating. We will use this comparison to fix the four free pa- 
rameters of our satellite orbit model, /orb, /A,h, /A,d and eh- 
As described in Appendix A, /orb controls the timescale on 
which tidally stripped mass is lost from the satellite, while 
/A,h and /A,d are the factors that appear in the Coulomb 
logarithms. Ah and Ad, for the dynamical friction force due 
to the halo and disk respectively. The parameter eh controls 
the strength of gravitational shock-heating and is defined in 
Benson et al. (2002a). 

Using our satellite orbit model, each orbit is integrated 
for 4 Gyr. Figures 2 and 3 show the orbital position and ve- 
locity and the remaining bound mass and orbital energy of 



the satellite for models G2S2 and G1S3 respectively, with 
our N-body results shown as open circles. Figure 3 also 
shows the energy deposited in the disk in model G1S3. This 
is given by the vertical kinetic energy of the simulated disk 
minus the vertical kinetic energy of the disk in model GISO 
which contains no satellite. The subtraction removes both 
the initial energy of the disk, and the energy gained by two- 
body relaxation during the simulation. We indicate at the 
top of each figure the label of the satellite model, the ini- 
tial inclination of the orbit with respect to the galaxy disk 
{9i), the initial circularity (e,j; the angular momentum of the 
satellite divided by the angular momentum of a circular or- 
bit with the same energy) and the initial apocentric distance 
of the orbit (ra). Where they are available, we show the re- 
sults of VW as triangles. Note that in the simulation of VW 
the satellite loses mass more rapidly, due to the incorrect 
density profile used. For comparison, we show, as dashed 
lines, the orbital radius and remaining bound mass derived 
from the analytical calculations of Taylor & Babul (2001) 
for the same model. 

The results in Figs. 2 and 3 are for the parameter com- 
bination (/orb, /A,h, /A,d, £h) = (2.5,1.5,3.0,1.0). The valucs 
of /orb and eh are fixed by matching the mass loss rates 
found in the simulations with no disk component. The value 
of /a.Ii, which controls the strength of the dynamical friction 
force due to the halo, is fixed by matching the rate of decay 
of the orbital radius in models with no disk (so that the or- 
bital decay is caused entirely by the halo plus bulge system). 
Finally, /A,d is fixed by matching the rate of orbital decay in 
the models which include a disk. In these models, the disk 
is the dominant source of dynamical friction throughout a 
substantial fraction of the orbital evolution. 

The parameter values that we have selected produce 
the best agreement with the set of fifteen models that were 
simulated. Generally, we find quite good agreement with the 
numerical results, comparable to that achieved by Taylor & 
Babul (2001)^. Our model uses more general expressions 
for Ah and Ad than that of Taylor & Babul (2001). If we 
treat those numbers as free parameters (instead of /A_h and 
/A,d) we are able to achieve even better agreement with the 
numerical simulations. However, our approach has the ad- 
vantage that Ad and Ah scale in a physically reasonable way 
when we apply our model to very different satellite/host sys- 
tems. In any case, orbital positions and velocities are typi- 
cally matched accurately until the final merging of the satel- 
lite (where it becomes difficult to determine these quantities 
precisely in the N-body simulations). The satellite mass as 
a function of time is typically matched to within about 30- 
40%. Table 3 lists several quantities - the final change in 
the disk energy and the time at which the satellite reaches 
50% and 10% of its original mass - from both analytic and 
N-body calculations for comparison. 

It is worth noting that the free parameters of our orbit 
model are set without reference to the disk heating rate seen 

■t It should be noted that Taylor &; Babul (2001) were attempting 
to match the simulations of VW, rather than our simulations, so 
that one should not expect exact agreement of their results with 
ours. 
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Model G2S2 (SI, 6*; = 0°, ej = 0.55, r-a = 55kpc, no disk) 




1 10 100 1000 10< 1 10 100 1000 

l/Gyr l/Gyr 



Figure 2. Evolution of the satellite and its orbit in the diskless model G2S2. We compare the results from our analytical model (solid 
lines) with the analytical model of Taylor & Babul (2001) (dashed lines) and with our N-body simulation (circles). Top left-hand panel: 
The orbital position and radius of the satellite as a function of time. Top right-hand panel: The orbital velocity of the satellite and its 
components as a function of time. Lower left-hand panel: The remaining bound mass of the satellite as a function of time. The dotted line 
shows the mass of the satellite if mass loss beyond the tidal radius is assumed to occur instantaneously (i.e. /orb = 0). Lower right-hand 
panel: The change in the specific orbital energy of the satellite as a function of time. 



in the numerical simulations. Thus, the heating rates we 
predict are entirely specified by other considerations. The 
lower right-hand panel in Fig. 3 shows the change in disk 
vertical kinetic energy from our analytic model calculated 
as described in §2.2.2 and from the N-body simulation. We 
find that our analytic model reproduces the final disk energy 
in the numerical simulations to better than a factor of two 
in ten out of the fifteen simulations (see Table 3) but, in 
extreme cases, the difference can be a factor of three or 
more. Of the five models which do not agree to within a 
factor of two, one (G1S7) has a prograde satellite orbit in 
the disk plane {di = 0°), two (G1S4 and GlSll) have polar 
orbits {Oi = 90°), and two (G1S8 and GISIO) are on inclined 
prograde orbits {6i = 45°). 



For all five of the most discrepant cases, the analytical 
calculation predicts less heating than the N-body simula- 
tion. The largest disagreement occurs for model G1S7 which 
has a prograde orbit in the disk plane. Here, the analytical 
determination overestimates the dynamical friction force in 
the disk as measured in the N-body simulation. The satel- 
lite then becomes trapped in an orbit rotating with the disk 
and there is no further energy transfer to the disk, resulting 
in and underestimate of the heating in the analytic model 
by a factor of 8. For the polar orbits (G1S4 and GlSll), 
mass loss in the analytic model is too rapid and this again 
reduces the heating rate compared to the N-body calcula- 
tion. These two models underpredict the N-body heating by 
a factor of approximately 3. For the inclined orbits (G1S8 
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Model G1S3 (SI, 6; = 45°, ej = 0.55, = 55kpc) 




t/Gyr 

Figure 3. 
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Figure 3. (cont.) Properties of tlie orbiting satellite and host halo galaxy disk in model G1S3. We compare the results from our analytic 
calculations (solid linos) and those of Taylor & Babul (2001) (dashed lines), with those from our N-body simulation (circles) and those 
of Velazquez & White (1999) (triangles). Top left-hand panel: The orbital position and radius of the satellite as a function of time. Top 
right-hand panel: The orbital velocity and speed of the satellite as a function of time. Centre left-hand panel: The remaining bound mass 
of the satellite as a function of time. The dotted line shows the mass of the satellite if mass loss beyond the tidal radius is assumed to 
occur instantaneously. Centre right-hand panel: The change in specific orbital energy of the satellite with time. Lower right-hand panel: 
The vertical kinetic energy of the central galaxy disk. Filled symbols show the energy measured in the original coordinate frame of the 
disk, whereas the open symbols show the energy measured in a frame that coincides with the principal axes of the inertia tensor of the 
disk at each epoch. The dotted line shows the result obtained if the energies of the disk in each direction (R, rp, z) are assumed to reach 
equipartition. 
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Figure 4. A comparison between the analytic and N-body re- 
sults for the change in the vertical component of the disk kinetic 
energy, ATz. The dashed line is the locus of perfect agreement, 
with the two dotted lines indicating factor of two discrepancies. 
The symbols, one for each of the fifteen simulations, GlSl to 
G1S15, indicate through their orientation, shape and shading the 
orbital inclination, the orbital circularity and the satellite model 
respectively. 

and GISIO), it is possible that the disk is no longer well 
described by a single inclination (for example, it may have 
become warped), leading to an overestimate of the energy 
in the N-body simulations. (With the number of particles 
employed in our simulated disks, the inclination of an un- 
warped disk can be determined to very high precision, so 
there is very little inaccuracy in the determination of the 
disk energy). 

Figure 4 compares the N-body and analytic results for 
the change, ATz, in disk kinetic energy. The dashed line is 
the locus of perfect agreement between the two calculations, 
and the two dotted lines indicate a factor of two discrepancy. 
The symbols, one for each of the fifteen simulations, GlSl to 
G1S15, indicate through their orientation, shape and shad- 
ing the orbital inclination, the orbital circularity and the 
satellite model respectively (as indicated by the key in the 
figure) . 

Several of the N-body simulations show evidence of bar 



formation in the central regions of the disk. This is particu- 
larly evident when the satellite is on a prograde orbit in the 
disk plane. Bars may be expected to enhance the transfer 
of energy to the disk, and may be part of the reason why 
the analytic model (which does not allow for bar formation) 
substantially underpredicts the amount of heating in some 
cases (e.g. G1S5 and G1S7, of which the latter shows a par- 
ticularly strong bar in the N-body simulation). 

The efficiency of vertical heating, e^, is an important 
component of our calculations. If we did not include this 
efficiency factor, the predicted heating rates would be up to 
4 times higher (depending on the orbit — the effect is largest 
for near-circular prograde orbits in the disk plane and polar 
orbits), with a factor of 3 being typical. 

The inclusion of the 8 dependence in the expression 
for Ad (see Appendix B2.2) tends to reduce the heating 
rate slightly. The effect is small for most orbits, but it is 
of greater importance for orbits in the disk plane, helping to 
improve the agreement with the simulations in those cases. 
The use of an anisotropic disk velocity dispersion in the dy- 
namical friction force generally has an even smaller effect, 
typically increasing the disk heating rate by a few percent 
(although in some cases the rate is decreased by an equally 
small amount). Prograde orbits in the disk plane are, once 
again, most strongly affected, with heating rates reduced by 
20-40% 

The galaxy in the N-body model contains a bulge of 
mass 1/3 that of the disk. VW also performed simulations 
with bulges of mass 1/5 and 2/3 that of the disk to examine 
the influence of the bulge on the heating rate, finding that 
a more massive bulge reduced the amount of disk heating. 
Our analytical model typically reproduces this trend, with 
approximately the same strength. 

To summarize, we are able to reproduce the rates of 
disk heating seen in numerical simulations for the majority 
of the cases considered. Where the analytic approach "fails" 
(we say "fails" , since the N-body techniques have their own 
inadequacies and do not necessarily represent the correct 
solution), it underestimates the heating by a factor of 3 on 
average. In many of the discrepant cases, the incorrect heat- 
ing rate is a consequence of an incorrect estimate of the disk 
dynamical friction force or tidal mass loss rate, but in some 
of the other discrepant cases, the reason is less obvious. It is 
worth emphasizing that our analytic calculation reproduces 
several important trends observed in the N-body heating 
rates. For example: 
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• heating is greatly suppressed for satellites on polar or- 
bits; 

• differences between heating rates for prograde and ret- 
rograde orbits (which are not always in the same sense, de- 
pending on the satellite type) are reproduced; 

• differences due to the concentration or initial mass of 
the satellite are clearly reproduced. 

The trend of increased heating for more circular orbits, seen 
in the N-body simulations, is not reproduced, however. 

While it is clear that the analytic model does not match 
the N-body heating rates perfectly, in the majority of cases, 
the differences are compatible with the accuracy of the sim- 
ulations themselves, as judged by the convergence tests. We 
conclude that, in general, the analytic model provides a rea- 
sonable approximation to the simulation results. 



4 RESULTS 

4.1 Scale-Height Distribution for Disk Galaxies 

Having demonstrated that our model can be used to calcu- 
late disk heating rates with reasonable accuracy, we now pro- 
ceed to apply these calculations to galaxy formation in a cos- 
mological setting. Specifically, we implement this model of 
disk heating in the GALFORM semi-analytic model of galaxy 
formation described by Cole et al. (2000) and Benson et 
al. (2002a), based on a standard ACDM cosmology with 

Qo = 0.3 and Aq = 0.7§. This model follows the growth of 
galactic disks in a merging hierarchy of dark matter halos. 
At each time the model predicts the mass and radial size 
of the galactic disk forming at the centre of each halo. It 
also gives the rate at which subhalos are merging into each 
halo, which we take as input for our calculations of satellite 
evolution and disk heating. We assume that only direct pro- 
genitors of the halo cause heating (i.e. subhalos can heat the 
disk, but sub-subhalos are not considered). This is to avoid 
double-counting of mass. We will consider briefly below the 
effect of allowing all progenitors to heat disks. 

Using this model, we generate a representative sam- 
ple of galaxies living in dark matter halos spanning a wide 
range of masses. For each galaxy, the model computes the 
usual properties predicted by this type of modeling (masses, 
luminosities, etc. — see Cole et al. 2000), and now also 
the vertical scale-height of the galactic disk. Figures 5 and 
6 show the resulting distribution of disk scale-heights, ex- 
pressed in units of the disk radial scale-length, for galaxies 
with Mb — 51og/i < —19.5 (approximately L, and brighter 
galaxies) and —19.5 < Mb — 51og/i < —17.0 respectively. We 
include only spiral galaxies (which we define as any galaxy 
with a bulge-to-total ratio measured in dust-extinguished 
B-band light less than 0.4). 

We remind the reader that we define the dimension- 
less scale-height, h = Hd/Rd, in terms of the thickness pa- 
rameter Hd in the sech^ vertical density law and the ra- 

§ Benson ot al. (2002c) describe small changes in the parameters 
of this model, relative to those of Benson et al. (2002a), which 
we retain here. 



dial exponential scale- length, Rd. The disk thickness can 
be equivalently defined as Hd ~ Ed/(2po), where Ed is 
the disk surface density and po the density at the mid- 
plane. However, many authors prefer to use the exponen- 
tial scale-height as the measure of disk thickness. Since 
sech^{z/Hd) oc exp{—2z/Hd) for z » Hd, the exponen- 
tial scale-height that would be measured for our assumed 
vertical profile is Hd,cxp = Hd/2. 

The left and right-hand panels in Figures 5 and 6 show 
the scale-height distributions for the global and local heating 
assumptions, from which we see that the results are not very 
sensitive to this choice. The figures also show the sensitivity 
of the results to two other parameters, one numerical and 
the other physical. 

The numerical parameter characterizes the galactocen- 
tric radius at which the satellite is assumed to merge with 
the main galaxy and stop heating the disk. In Cole et al. 
(2000) and in Paper I, we assumed that two galaxies merge 
at the time when the separation of their centres, 7?, equals 
the sum of their half-mass radii, Ri^^^ -\- Ri^^^^^. However, 
once tidal stripping is taken into account, it would seem 
reasonable to allow the satellite to sink down to 7? = 
and continue heating the disk while it does so. However, 
for numerical reasons we cannot integrate the satellite or- 
bits down to i? = 0. We have therefore calculated the 
disk heating when the satellite orbit is followed down to 
R = /heat(-Ri,3.t + ^ihost) for /heat = 1, 3 and f. We 
find that the distribution of scale-heights has converged for 
/heat = |, and use this as our standard value in what fol- 
lows. We show in Fig. 5 results for /heat = 1 (dot-dashed 
lines) and /heat = | (solid lines) in both cases with disk 
heating by molecular clouds also included. The differences 
in the scale-height distributions are fairly small (they are 
somewhat more significant if we do not include disk heating 
by molecular clouds). 

The physical parameter concerns the heating of the 
disk that results from scattering of stars by giant molec- 
ular clouds, computed using eqn. (9). Our standard calcula- 
tion (solid lines in Figs. 5 and 6) includes heating by clouds 
with the parameters described in §2.2.5. The figures also 
show the results when no clouds are present (dotted lines) 
and when the masses of individual clouds and the fraction of 
gas in clouds are both doubled (dashed lines). Removing the 
clouds entirely results in a tail to very low h in the height dis- 
tribution. These galaxies experienced very little heating by 
substructures, and so their thickness is almost entirely due 
to heating by molecular cloud scattering. The peak of the 
distribution is little changed, but the median scale-height is 
significantly reduced (see Figure 7). Doubling the cloud mass 
(dashed lines) results in a shift towards somewhat thicker 
disks without changing the shape of the distribution. 

Global and local heating are found to produce rather 
similar distributions of scale-heights. Note that we have 
shown the results for local heating at the disk half-mass 
radius. Our local heating model, in fact, predicts a trend of 
increasing scale-height with disk radius; we defer a detailed 
study of this radial variation to a future paper. 

It should be noted that the tails of the distributions ex- 
tend to h > 1, which is clearly unphysical. Our analytical 
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Figure 5. Normalized disk scale-height distributions for spiral galaxies with Mb — 51og/i < —19.5. Left and right-hand panels show 
results for global and local heating respectively. Solid lines show results from our full calculation, including heating from substructure 
and from scattering by molecular clouds. The dotted line corresponds to ignoring the molecular cloud heating, while the dashed line 
corresponds to increasing the masses of individual clouds and the total mass in clouds by a factor of two over our standard values. These 
results correspond to satellite orbits which are integrated until they reach a radius {^1^^^^ + ^iijogt)/8. For comparison, the dot-dashed 
line shows the result when the integration is stopped when a radius [Ri^^^ + R^^^^g^) is reached. The vertical shaded strip shows the 
observationally allowed range for the scale-height of the Milky Way galaxy, discussed in §4.4. 



calculation is based on the thin disk approximation, h 1, 
and so breaks down when /i ~ 1. We interpret these ob- 
jects as disks that have been heated so much that they are 
no longer disks, and must instead resemble a spheroidal or 
irregular galaxy. For these galaxies, our calculations break 
down, but we can safely assume that they are no longer 
recognizable disk galaxies. 



4.2 Other Effects on the Scale-Height 
Distribution 

We now discuss tests of various potential systematic effects 
in our calculations. 

Merger Tree Resolution: Our calculations typically re- 
solve dark matter substructures with mass greater than 
5 X 10^ Mq in every merger tree. Thus, we ignore the 
heating due to lower mass halos. Increasing the resolu- 
tion of to 10^ Mq results in no significant increase in 
the amount of heating experienced by galaxies, indicating 
that our resolution is sufficient to estimate the total heat- 
ing rate. (Note that the heating produced by a satellite of 
mass M should scale approximately as , making it rela- 
tively easy to achieve convergence here provided the number 
of satellites, dA''/dlnM, varies with mass less steeply than 



M ^ at small mass. In fact, numerical simulations indicate 
dN/dlnM ~ M'^-'' for subhalos (Springel et al. 2001).) 

Ejfects of Sub-subhalos: In our standard calculation, 
sub-subhalos (i.e. halos which reside inside a larger halo 
which subsequently fell into a yet larger halo) do not con- 
tribute separately to the heating of disks. (Note that this is 
different from our treatment of galaxy mergers; the merging 
times of sub-subhalos are computed from their own prop- 
erties, not those of the subhalo in which they reside.) An 
alternative approach would be to treat sub-subhalos (and 
higher levels of the merging hierarchy) on an equal basis as 
subhalos. To avoid double-counting of mass in this case, we 
must remove the mass bound to sub-subhalos when deter- 
mining the mass of a subhalo. We do this by scaling down the 
density profile of the subhalo so as to remove this amount 
of mass before computing heating rates. 

If we adopt this approach, we find that the distribution 
of scale heights is shifted to slightly larger values. Heating 
by sub-subhalos, however, would only be important if these 
sub-subhalos survived after their host had been tidally de- 
stroyed. While it is unlikely that this would occur to any 
great extent, numerical simulations could, in principle, an- 
swer this question. 

Effect of Cosmological Model: Finally, we have repeated 
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Figure 6. Normalized disk scale-height distributions for spiral galaxies with —19.5 < Mb — 5 log ft < —17.0. Left and right-hand panels 
show results for global and local heating respectively. Solid lines show results from our full calculation, including heating from substructure 
and from scattering by molecular clouds. The dotted line corresponds to ignoring the molecular cloud heating, while the dashed line 
corresponds to increasing the masses of individual clouds and the total mass in clouds by a factor of two over our standard values. 



our calculations in an fio = 1 cosmology, using the rCDM 
parameter set used by Benson et al. (2000a), but including 
the effects of photoionization suppression^. This model is 
not as successful at matching the properties of z = galax- 
ies as our standard ACDM model. In particular, galaxies 
are somewhat too faint to match the observed luminosity 
function (by about 0.75 magnitudes in the B-band), forcing 
us to adopt an unphysical value of the mass-to-light ratio 
normalization parameter, T of 0.7. We find that the median 
scale-height of L, disk galaxies is slightly smaller in this cos- 
mology than in our ACDM model. At first sight, this seems 
surprising, since, as noted by TO, there is more infall of sub- 
structure at late times in an fio = 1 cosmology, which would 
result in a larger rate of heating at the present day. However, 
our model galaxies in this cosmology are younger than their 
ACDM counterparts (due to the later growth of structure 
and to the stronger feedback required in this model) , and so 
they have less time in which to be heated. These two effects 
counteract each other. 

The age of our Galactic disk has been estimated using 



^ Note that Benson et al. (2000a) adopted an artificially high 
merger rate in order to obtain a good match to the galaxy lumi- 
nosity functions. With our more detailed model of merging, we no 
longer have the freedom to adjust the merger rate in this way. We 
find that, in this cosmology, our revised merger model produces 
somewhat too few elliptical galaxies. 



studies of white dwarfs. For example, Fontaine, Brassard & 
Bergeron (2001) find an age of 11 Gyr for the Galactic 
disk. Since they assume a constant star formation rate, this 
implies a mean stellar age of 5.5 Gyr for the disk stars. Our 
rCDM disks typically have a mean stellar age of 4 Gyr, 
somewhat less than the true value. Consequently, our model 
underestimates the amount of heating experienced by disks 
in this cosmology, albeit only by a small factor. 

We can understand the similarity of the disk scale- 
heights in the two cosmologies in more detail by examin- 
ing the growth histories of the dark matter halos hosting 
L* spiral galaxies. In our model, halos of present-day mass 
2 X lQ^^h~^ Mq have, on average, assembled half of their 
mass by redshifts of 0.45 and 0.91 respectively in the rCDM 
and ACDM cosmologies. The mean stellar ages of disk 
galaxies — 4.0 and 5.5Gyr for rCDM and ACDM respec- 
tively — reflect this difference in halo assembly epoch. We 
find that the host halos on average accrete close to 25% 
of their total mass over these galaxy lifetimes in both cos- 
mologies. Therefore, the number of substructures infalling 
onto a galaxy over its lifetime is roughly the same in both 
cosmologies, consistent with their similar disk scale-height 
distributions. It should be kept in mind that the disks in our 
Slo = 1 model are somewhat unrealistic (e.g. they are too 
faint for a reasonable T and, more importantly, too young). 
An r^o = 1 model which produced realistic disks might pre- 
dict larger (or smaller) scale-heights. The important lesson 
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Figure 7. The median fractional scale-heights, h = H^/R^, of 
spiral galaxies as a function of absolute magnitude. The squares 
show the results for heating by substructures only, and the circles 
for heating by substructures and clouds together. In each case, the 
filled symbols arc for global heating, and the open symbols (offset 
slightly for clarity) for local heating. The errorbars indicate the 
10% and 90% intervals of the distribution of scale-heights. For 
clarity, the errorbars arc suppressed for the local heating case, 
but are similar to those for global heating. 



to derive from these considerations is that disk scale-heights 
depend on the details of galaxy formation as well as on the 
cosmological model. 



4.3 Scale-Heights as a Function of Luminosity 

In Fig. 7, we show the median value of ^ as a function of ab- 
solute magnitude for spiral galaxies in our standard model. 
The squares show the results for heating by substructure 
alone, and the circles for heating by substructure and clouds 
together. The median scale-height at all luminosities is much 
larger when heating by clouds is included. However, the scat- 
ter in scale-height at a given luminosity is extremely large 
for the case of heating by substructures only, reflecting the 
strongly stochastic nature of this process. Our calculations 
predict that brighter galaxies should host thinner disks than 
fainter galaxies (when measured in terms of the fractional 
disk thickness, h = H^/Rd)- This trend is apparent in cal- 
culations with and without molecular clouds, and reflects a 
similar trend in the fractional vertical energy. 



4.4 Comparison with the Milky Way galaxy 

It has been conventional to compare predictions for disk 
scale-heights with the observed value for the Milky Way 
galaxy. As a way of testing models against the real uni- 
verse, this comparison has significant drawbacks, since (i) 
the global parameters of the Milky Way (such as the disk 



radial scale-length, total luminosity and bulge-to-disk ra- 
tio) are, in fact, quite difficult to determine observationally 
and (ii) the models predict a distribution of scale-heights 
at a given luminosity, and this cannot be constrained well 
from a single measured point. Therefore, we will make only 
a brief comparison with the Milky Way here, before com- 
paring with the distribution of scale-heights measured for 
external galaxies. 

The vertical scale-height of the galactic disk in the So- 
lar neighbourhood has been measured from star counts. 
We use the recent determination by Mendez & Guzman 
(1998) which, for a sech^ {z / Hd) vertical profile, gives Hd ~ 
0.50±0.08kpc (corresponding to an exponential scale-height 
of 0.25 ± 0.04kpc), somewhat smaller than earlier determi- 
nations. The measurement of the radial exponential scale- 
length of the galactic disk has, in the past, been a matter 
of more disagreement. We use the models of the galactic 
mass distribution by Dehnen & Binney (1998), which imply 
Rd ~ 3.0 ± 0.4kpc. Combining these, we find the fractional 
scale- height of the Milky Way stellar disk, h = Hd/Rd = 
0.18 ± 0.05. This range in h is indicated as a shaded region 
in Figure 5, from which one can see that the scale-height of 
the Milky Way is entirely typical of L* disk galaxies in the 
model (with 35% of galaxies predicted to have h > 0.18). 
We have repeated the comparison using the same definition 
of "Milky Way-like" galaxies as in Benson et al. (2002b), 
namely a circular velocity at the disk half-mass radius be- 
tween 210-230km/s and a bulge-to-total ratio by mass be- 
tween 5-20%. We again find that the observed scale-height 
of the Milky Way lies well within the distribution of h pre- 
dicted by the model (with 80% of such galaxies predicted to 
have h > 0.18). 



4.5 Comparison with the Observed Scale-Height 
Distribution for Other Galaxies 

The best way to test models of disk heating is by comparing 
with the observed distribution of scale-heights for external 
galaxies. This distribution has recently been measured in 
a complete sample of disk galaxies, for the first time, by 
Bizyaev & Mitronova (2002). They estimated the vertical 
and radial scale-lengths of a statistically complete sample 
of 60 edge-on galaxies using K-band photometry from the 
2MASS survey. 

We compare the scale-height distribution in Bizyaev & 
Mitronova (2002) sample with our model predictions in 
Figure 8. Since the selection criteria for the observational 
sample are somewhat complex, we weight model galaxies 
so as to match the distribution of absolute magnitudes in 
the observational sample (which peaks in the range — 19 < 
Mb — 51og/i < —18), and select only those galaxies with 
bulge-to-total luminosity ratios typical of the morphological 
types found in the observational sample (which are mostly 
Sc spirals). We see that the model provides quite a good 
match to the observed distribution, with both distributions 
peaking around h = 0.2. The only significant discrepancy is 
that the model predicts too many systems with large h ^0.4. 
However, it is not clear that such thick galaxies would be rec- 
ognized as disk galaxies. The conclusions that can be drawn 
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Figure 8. The normalized distribution of scale-heiglits, h, in the 
observational sample of Bizyaev & Mitronova (2002) compared to 
the prediction of our model. Error bars on the observational dat- 
apoints indicate Poisson errors. The model predictions are shown 
by solid and dashed lines for global and local heating respectively. 
The model galaxies have been weighted to match the distribution 
of absolute magnitudes and morphological types found in the ob- 
servational sample. 



at present are limited by the relatively small size of the cur- 
rent observational sample. However, this situation should 
soon improve with the availability of data from large CCD- 
based sky surveys, which will allow much more thorough 
tests of the theoretical predictions. 



5 DISCUSSION 

We have developed a model to calculate the rate of heating 
of galactic disks by substructures orbiting in their halos. 
To calibrate the model, we performed N-body simulations 
of disk heating which we tested for convergence. We find 
that the analytical model reproduces the heating rates in 
the N-body simulations to within a factor 3 in most cases. 
One could perhaps improve the accuracy of the analytical 
model by treating the satellite-disk interaction in terms of a 
sum of interactions with resonances in the disk (e.g. Donner 

6 Sundehus 1993; Weinberg & Katz 2002). It is unclear, 
however, whether such a calculation would ever be worth 
performing semi-analytically, i.e. whether its computational 
cost would be any less than that of a full N-body simulation. 
Nevertheless, it is clear from the calculations presented here 
that N-body estimates of disk heating rates have their own 
problems (e.g. very large numbers of particles are required 
in the disk to determine the heating rate accurately) , and so 
it may yet prove worthwhile to pursue analytical estimates 
of disk heating. 

We find that for galaxy formation in the standard 



ACDM cosmology, heating by substructure alone produces 
a distribution of disk scale-heights which is very broad and 
skewed to low values, with median fractional scale-height, 
h = Hd/Rd, around 0.05 for L* spiral galaxies. The width 
of the distribution reflects the stochastic nature of the heat- 
ing process, which is, in turn, related to the distribution 
of orbital parameters of the satellites. Including the addi- 
tional heating generated by stars scattering from gas clouds 
in the disk increases the median value of h significantly, to 
around 0.2. The distribution is considerably less broad once 
the contribution from gas cloud heating is included. Heat- 
ing by clouds is treated as a deterministic process here, with 
variations in the amount of heating for a given type of galaxy 
reflecting the distribution of ages of galactic disks. The frac- 
tional scale-height for the Milky Way galaxy, estimated ob- 
servationally to be around 0.2, is then entirely consistent 
with our model expectations for a typical L, spiral galaxy. 
We find that the predicted distribution of scale-heights for 
slightly sub-L* spiral galaxies agrees remarkably well with a 
recent observational determination by Bizyaev & Mitronova 
(2002) based on data from the 2MASS survey. 

It is intriguing that for the i, galaxies considered here, 
satellites and gas clouds give rise to comparable amounts 
of disk heating. A simple order of magnitude estimate of 
the scale-heights produced by these two processes illustrates 
why this is so. Using the expressions given in this paper, we 
find (for parameter values typical of Milky Way-like galax- 
ies) that the fractional scale-height generated by scattering 
from giant molecular clouds is 



h 



7.2 X 10" 



3 


r / 1 


1/2 


■ MJMd ' 


1/2 


MnAl 


1/2 




L0.025J 




3 X 10-5 




[ 3 J 





V 


1/2 


«s(/?)' 


3/2 


■^s(/3)' 


1/2 


t 


90 Gyr-\ 




0.7 




0.15 




Gyr 



1 1/2 
(11) 



while that generated by dark matter substructures is 
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In the first equation, / is the fraction of the total disk 
mass in the form of giant molecular clouds. In the second 
equation, 14 is the typical orbital velocity of satellites, ro^b 
their typical orbital radius, /mass is the fraction of the to- 
tal halo mass in the form of substructures, /max is the mass 
of the largest substructure in units of the total halo mass 
and we have assumed a distribution of substructure masses 
AN I AM oc M-^-\ In both cases t is the time for which heat- 
ing has occurred. To derive the second expression, we have 
assumed that substructures heat the disk only over a frac- 
tion of their orbit approximately equal to Hd/roA- Taking 

t ~ 10 Gyr, these estimates imply h ~ 0.1 1 for both 

heating mechanisms, confirming the coincidence that the 
two contribute approximately equally to the scale-height of 
Milky Way-like disks (given the crude approximations made 
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above and the fact that we have ignored the stochastic na- 
ture of heating by satelhtes). However, these two expressions 
have different dependencies upon the properties of the galax- 
ies in question. Thus, we should not expect the two to make 
equal contributions to the scale-height of galaxies dissimilar 
to the Milky Way. This may be seen in Fig. 7, where it is 
clear that the heating by substructures is relatively less im- 
portant for lower luminosity galaxies. In conclusion, the fact 
that the two heating mechanisms make similar contributions 
to the scale-heights of Milky Way-like galaxies appears to be 
coincidental. 

It is interesting to compare our conclusions with those of 
TO, who found that the Milky Way disk could have accreted 
only up to 5% of its mass within the Solar circle within 
the past 5 Gyr without becoming too thick. Our calculations 
show that the Milky Way halo in fact accreted around 25% 
of its mass (i.e. the total dark mass of the halo) during 
this time. This is approximately 100 times more than the 
TO limit. Disks in our model are able to remain fairly thin 
despite this substantial accretion for two reasons. Firstly, 
many of the accreted subhalos have orbits which do not take 
them close to the central galaxy disk, and so they contribute 
almost nothing to the heating of the disk. Figure 9 shows 
the amount of energy transferred to the disk of a Milky 
Way-like galaxy by individual dark matter satellites as a 
function of their mass. At a fixed mass, the distribution of 
heating energies has a bimodal distribution. Virtually all 
of the heating energy is supplied by the satellites in the 
upper branch which are those whose orbits take them close 
to the central galaxy disk. These satellites are "trapped" by 
dynamical friction and damage the disk during an extended 
period; distant satellites, on the other hand, have a negligible 
effect. The satellites that cause most of the heating amount 
to only 6% by mass. Thus, of the 5 x lO^M© infalling, only 
around 3 x 10^" Mq contribute to heating the disk. This is 
still 6 times larger than the TO limit (0.5 x lO^^M©). 

The second difference with TO is that we find that tidal 
mass loss in subhalos substantially reduces the amount of 
heating experienced by the disk. This is in disagreement 
with TO, who found that tidal mass loss reduced disk scale- 
heights by at most a factor of two. If we do not allow satel- 
lites to lose mass, the peak of the scale-height distribution 
is shifted to a value of h which is approximately ten times 
larger than our standard result. This, of course, reflects the 
different density profiles that we assign to both the host and 
satellite halos (and which are significantly more extended 
than the objects considered by TO), and the associated in- 
crease in the dynamical friction timescale in our model. In 
conclusion, our results are in partial agreement with those 
of TO — halos in our model accrete much more mass in the 
past 5 Gyr than the TO limit, but little of this mass ever 
contributes to heating the disk. 

There is clearly a need for further study of the heating 
of galactic disks. In particular, the importance of heating by 
satellite-triggered bars and the extent to which heating is lo- 
cal or global are important, yet poorly understood aspects of 
the problem. We believe that analytical modeling of the type 
developed in this paper provides a powerful means by which 
to estimate the degree of heating by substructures and could 




Figure 9. The energy contributed to disk heating by satellites 
as a function of their mass in shown in the lower panel. Points 
show the results for a large sample of satellites, for which we plot 
the heating energy supplied {AE^ disk) against the mass of the 
satellite (expressed in units of the mass of the host system halo). 
At fixed mass, the distribution shows a bimodal form, the dashed 
line indicates the approximate division between the two peaks of 
the distribution. Crosses indicate the mean heating energy per 
satellite at each mass, and the dot-dashed line shows an approxi- 
mate fit to these points. In the upper panel, we show the fraction 
of points at each mass which lie above the dashed line in the lower 
panel. 



easily incorporate any improvements in our understanding 
of the physics of the process. Its particular strengths are the 
ability to resolve fully all substructures contributing to the 
heating and to compute many realizations of the heating 
process rapidly, thus allowing the full distribution of scale- 
heights to be determined. These features have allowed us to 
present predicted distributions of galaxy scale-heights which 
will be tested by forthcoming observational data. 

In conclusion, the observed thickness of the Milky Way's 
stellar disk seems to be entirely consistent with the amount 
of substructure in galactic halos expected in a cold dark mat- 
ter universe. Stars scattering from giant molecular clouds 
and substructures passing through or near the disk pro- 
duce similar amounts of heating. Distinguishing these two 
contributions observationally might be possible by means of 
the stellar age-velocity dispersion relation in the Milky Way 
disk. An important extension of this work will therefore be 
to examine model predictions for heating as a function of 
time within individual galaxies. The lowest values of h for 
disk galaxies are set by the heating due to star-cloud inter- 
actions, while the highest values are set by the heating due 
to substructures. Thus, precise measurements of the disk 
scale-height distribution can potentially constrain these two 
processes. 
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APPENDIX A: IMPROVEMENTS IN THE 
SATELLITE EVOLUTION MODEL 

We detail here the changes in the model of satellite evolution 
presented in Paper I. 

i) As before, that mass of a satellite which has become 
unbound due to tidal forces is lost gradually over a time com- 
parable to the orbital period. The fraction of the unbound 
mass lost in a small timestep of duration St was chosen to 
be proportional to St/torh, where torb is an estimate of the 
orbital timescale. In Paper I, we chose torb = 2'k/u}, where 
ul> is the instantaneous angular velocity of the satellite. In 
the present work, we instead take forb = /orb 27r/.^LjpcriWapo , 
where and Wapo are the angular velocity of the satel- 
lite at its most recent pericentric and apocentric passages, 
and /orb is an adjustable parameter which we expect to be 
of order unity. (Prior to the first pericentric passage we re- 
vert to our previous definition of mass loss timescale — this 
makes little difference to our results as typically very little 
mass loss occurs prior to this time.) The advantage of this 
choice is that it produces smoother mass loss histories (as is 
shown in §4). Furthermore, when considering cosmological 
distributions of satellites, we occasionally find orbits which 
are near-radial. The above definition then prevents the mass 
loss rate from becoming arbitrarily small. 

ii) In Chandrasekhar's formula for the dynamical fric- 
tion force, Taylor & Babul (2001) adopted fixed Coulomb 
logarithms of In A = 2.4 for the dynamical friction force due 
to the combined halo/bulge system and InA = 0.5 for the 
force due to the disk. They found that these values resulted 
in the best match to the results of the numerical simula- 
tions of Velazquez & White (1999; hereafter VW), and we 
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adopted the same values in Paper I. Since we will be inter- 
ested here in a wide range of satellite and host halo masses, 
we adopt more general definitions. For the halo and bulge 
systems we take Ah — /A,h''(«sat + ''"iD)/GMsat, where r is 
the orbital radius of the satellite, Meat its mass, iigat the 
orbital velocity of the satellite, ctid the one dimensional ve- 
locity dispersion of the halo at radius r, and /A,h is a pa- 
rameter. Since Ah < 1 is possible with this definition, we 
replace the usual In Ah term in the expression for the dy- 
namical friction force (eqn. 20 of Paper I) with i ln(l -|- Ah), 
the correct form for small Ah (Binney & Tremaine 1987). 
We also account for the finite size of the satellite as described 
in Appendix B2. For the disk we must account for the dif- 
fering scale-lengths in the radial and vertical directions. A 
suitable expression for the Coulomb logarithm is derived in 
Appendix B2.2, and depends on the disk scale-length and 
velocity dispersions, the velocity of the satellite relative to 
the disk, the angle this velocity makes with the disk plane, 
and on a parameter, /A,d which plays a similar role to /A,h. 
These forms are used throughout our calculations. 

iii) The disk is now treated as having an anisotropic 
velocity dispersion (ctr, cr^, o"z) in the radial, azimuthal and 
vertical directions, and this anisotropy is included in the 
calculation of the dynamical friction force due to the disk 
(see Appendix B2.1). We adopt essentially the same model 
for the disk velocity dispersion components as VW. For 
the radial velocity dispersions, we set cr^ oc exp{—R/Rd) 

(Lewis & Freeman 1989)11 , where Rd is the disk radial scale- 
length, and fix the normalization by assuming the disk to 
have a Toomre Q-parameter of 1.5 at its half-mass radius, 
which results in Q ~ 1.5 at the Solar radius in a Milky 
Way-like galaxy disk (VW). The azimuthal velocity disper- 
sion is then determined using the epicyclic approximation, 
CT^ = a^K^/4Q^ (where k is the epicyclic frequency and Q, 
the orbital frequency of the disk). The vertical velocity dis- 
persion at each radius is calculated from the vertical scale- 
height Hd, assumed constant with radius, using the expres- 
sions in §2.2.2 (the vertical scale-height in turn is related to 
the disk vertical energy). In the analytical disk-heating cal- 
culation, the radial and azimuthal velocity dispersions are 
kept fixed in time, but the vertical velocity dispersion evolves 
with the disk vertical energy. 

iv) When computing the dynamical friction force due to 
the disk, we smooth the disk density to account for the finite 
size of the satellite halo as did Taylor & Babul (2001). We 
smooth on a scale equal to the current radius of the satellite 
after tidal limitation and gravitational shock-heating. 

v) As the disk scale-height will increase as a function of 
time due to disk heating, we allow for a variable disk scale- 
height in our satellite orbit calculations. This affects both 
the dynamical friction force due to the disk and also the 
gravitational forces exerted by the disk. 

vi) Heating by gravitational shocks causes shells of ma- 
terial within a satellite to expand before they become com- 
pletely unbound. Previously, this effect was included in the 



II Note that VW contains an error in this equation, although the 
text of that paper is correct. 




Figure Bl. The geometry of a scattering event contributing to 
the dynamical friction force on mass M (the satellite) due to a 
particle of mass m (a background particle). Here, Vm and are 
the velocities of M and m respectively, while Vq is the relative ve- 
locity of the two and b is the impact parameter for this scattering 
event. 

calculation of the tidal mass loss, but not in the calcula- 
tion of the final internal structure. We now calculate the 
evolution of the internal density and circular velocity pro- 
file assuming that the radii of shells of dark matter scale 
in inverse proportion to their energy. We have repeated the 
comparison we performed in Paper I of the distribution of 
peak internal circular velocities of satellite halos predicted 
by the semi-analytical model with the results of cosmolog- 
ical N-body simulations. We find that the same choice of 
initial satellite orbital parameters as in Paper I still gives 
the best match to the N-body simulations. 



APPENDIX B: DYNAMICAL FRICTION 
FORMULAE 

In this appendix we derive several formulae related to dy- 
namical friction which are employed in this work. For com- 
pleteness, in §B1 and §B2 we derive several well-known re- 
lations relevant to dynamical friction. A more complete dis- 
cussion of these results can be found in Binney & Tremaine 
(1987) for example. We consider a mass M moving through 
an infinite and homogeneous sea of particles of mass m{<^ 
AI), number density n and density p = mn. 

Bl Single Scattering Events 

For a single scattering event we take the results of Binney 
& Tremaine (1987; page 422). The scattering geometry is 
illustrated in Fig. Bl. The changes in M's velocity parallel 
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and perpendicular to the initial relative velocity vector of 
the m and M, Vo, are: 



AV| 
and 
AV, 



2m 
17 



1 + 



G2A/2 



Vo, 



2mVy 
GAP 



1 + 



(Bl) 



(B2) 



where b is the impact parameter and we have assumed the 
background particles to be much less massive than the object 
for which the force is being calculated. 

B2 Rate of Change of Velocity 

We now envisage a sea of particles m with a distribution 
of velocities given by /(Vm)- The contribution to the rate 
of change of velocity in the parallel direction from particles 
with velocity Vm is simply 



dV 



/(Vm) 



27r6nVoAV||d6, 



(B3) 



dt 

This gives, 

dVii , , Vn 

where A — femaxVo^/GM. If M has a finite extent (corre- 
sponding to replacing the lower integration limit of with 



(B4) 



by 

Acff 



1 + 



1 + [&min/femax]2A2 



- 1 



1/2 



(B5) 



Throughout this work, we take bmin equal to half the current 
tidal radius of the satellite. 

Clearly, the net change in the velocity of M perpendic- 
ular to Vo is zero by symmetry. Thus, the net rate of change 
of velocity of M is 

'^^^ 27rln(l+A2)pG^M / /(Vm) /J'" ' J^^^ dVm.(B6) 

Vm) 



dt 



(V„ 



The integral in the above equation has an identical form to 
integrals used to find the gravitational force at position 2:0 
due to a density distribution, if we identify /(Vm) = Gp(x), 
Vm = X and Vm = xq. 

Thus, the power extracted from the body through dy- 
namical friction is given by, 

dVM 



dt 



= 27rln(l-f A^)pG^A/^ 



, , (Vm - Vm) 3,^ 
'"^Vm-VM)3 " 



(B7) 



B2.1 Application to an Arbitrary Velocity Ellipsoid 

Binney (1977) derives an expression for the dynamical fric- 
tion force due to a system of particles with uniform density 



and Gaussian velocity distribution with dispersion a± in 
one direction and (T|| in the other two directions. Binney's 
equation (A4) is trivially generalized to the case where the 
velocity dispersions differ in all three directions. Combining 
this with his equation (A3) we find the following expression 
for the dynamical friction force: 



Fdf = 



2n HI + A')pGHi'^^^^^ 
X (Baf ReR, + B^v^e^ + B^v^e^) 



(B8) 



where p is the background density, M the mass of the or- 
biting object, {vn,v<t>,Vz) is the relative velocity vector of 
object and background particles (in cylindrical polar coor- 
dinates since we will apply this expression to a galaxy disk) , 
Br, 60, Gz are the basis vectors of the cylindrical polar co- 
ordinate system. The coefficients B axe given by. 



Br = 



dq 



[(l + g)3(l-e2-f g)(l-ei+g)]i/2 



X exp 



1 

~2 

2 /„2 



2 / 2 

(1 + g) ' (1 



+ 



^1 + 1) 



2 I A 



(1 



+ g) 



(B9) 



Bj, — 



dq 



[(l + g)(l-e2+ 5)3(1 _e2+g)]i/2 



X exp 



1 
2 

2 / 2 



+ 



(1 + g) (1 



e| + g) 



+ 



z / z 
«z/o-R 



(1 - ei + g) 



(BIO) 



Bz 



dq 



X exp 



[(H-g)(l-e2 +g)(l-ei+q)3]i/2 



1 
2 

2 /_2 



(1 + g) 



(1- 



^1 + 1) 



+ 



z I z 
"z/c^R 



(1 



q) 



(Bll) 



where 1 



4 = '^^/'^R and 1 - 



B2.2 Effective Coulomb Logarithm for the Disk 

In calculating the dynamical friction force due to the disk 
we require the Coulomb logarithm, | ln(l -I- A^), where A 
is normally defined as A = bmax^o^/GM, where Vo is the 
typical relative velocity of the satellite and stars in the disk. 
We adopt Vq = V,li + (o-| + + a^)/3, where Koi is the 
relative velocity of the satellite and the bulk disk motion, 
and ctr, (j0 and are the three components of the disk 
velocity dispersion. 

When computing the dynamical friction force we sum 
the contributions from all particles with impact parameter 
b by integrating around an annulus of radius b normal to 
the relative velocity vector of the particles and the satel- 
lite. We define tp as the angle of a point on this annulus 
measured from a vector, p, which lies in the plane of the an- 
nulus and which is parallel to the disk plane (see Fig. B2). 
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Figure B2. The geometry used in calculating the Coulomb loga- 
rithm for the disk. Scatterings through an annulus normal to Vq 
(the relative velocity of satellite and disk stars) are considered, p 
is a vector lying in this annulus and parallel to the galaxy disk. 
Angles in the annulus, 5/1, are measured from p. Finally, 9 is the 
angle between Vq and the normal to the disk. 



For the disk, the value of 6max, the upper limit of integration 
in eqn. (B3), will vary as a function of tp. In the direction 
corresponding to -0 = (and -0 = tt) the disk density distri- 
bution has a characteristic length-scale of (the exponen- 
tial scale- length). This will therefore correspond (approxi- 
mately) to the largest impact parameter scatterings occur- 
ring in that direction. In perpendicular directions {tp = tt/2 
and 1/) — 3n/2) a more appropriate characteristic length is 
reff = 7id(cos^ + sin^ ^)^''^, where 6 is the angle between 
the satellite-disk relative velocity vector and the z-axis, and 
h the ratio of disk scale-height to scale-length. Thus, the 
effective A in direction -0 is 



/A.dfidfe'V? 

GM 



(B12) 



where h' = (cos^ V + [cos^ 6 + sin^ sin^ V')^''^- If ac- 
count for the finite size of the satellite then: 



(1 + A')cff 



1 + (bminA/i?d/i7A,d)2 



(B13) 



The effective Coulomb logarithm is found by averaging over 
all ip: 

^iln(l-HA2),ffJ) = ^ / Hl + A^sd^p. (B14) 



This integral is solved numerically. 



Figure B3. Geometry used in computing the rate of increase of 
velocity dispersion in direction n. Vectors b, bi and b2 lie in the 
plane of the annulus. Vector Vq is normal to the annulus and 
vector n lies in the plane of Vq and bi . 



B3 Rate of Increase of Scattered Particle Velocity 
Dispersion 

We now wish to determine the rate of increase of the one- 
dimensional velocity dispersion, measured in direction n, 
of the particles m due to dynamical friction scatterings. 
Since the centre of mass remains fixed during the scatter- 
ing, mAVin + A4"AVm ~ 0. Therefore, to find the change in 
velocity of m we multiply the equations (7-lOa) and (7-lOb) 
of Binney & Tremaine by —M/m. Writing these in a more 
convenient form: 



-2V0A7 



AV;„|| =-2Vo 



1 + A' 



1-f A" 



b'L 



(B15) 



(B16) 



for the components of velocity perpendicular and parallel to 
the relative velocity vector Vq as measured in the frame in 
which the centre of mass of M and m is at rest. 

Consider now the velocity of m in the frame in which 
the centre of mass of the central galaxy and its halo is at 
rest. The velocity changes are independent of frame so the 
final velocity of m in this frame is: 



^i?=V„, + ^AV;,|| + ^AK, 
Vo "6 



(BIT) 



We are interested in the velocities in some direction n. The 
initial and final velocities of m in this direction are: 



m,n 



(B18) 

V,n • n + AV,-n\\ cos 61^0 + AVm± cos St, (B19) 
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where 6vo and 9i, are the angles between n and Vq and fi and 
b respectively. The change in the component of the kinetic 
energy in direction n is therefore 



AEa = y {AVl^ I cos" ev„ + AVl^ cos" 
+2AVni||A'l4i± cos^Vq cos 9b 

+ 2Vn, • n [AVn-,\\ cosOvo + AVn, 



I } (B20) 



To sum over all particles m, we first integrate around an 
annulus of constant |b|. On this annulus, Ovq is constant 
and we can write the vector b = bi cos (f> + h2 sin 0, where 
bi • b2 = 0, |bi| = |b2| = b and is a parameter (see 
Fig. B3). We then note that 



cos 9bd4> — 



/'2ir 


bi • n 


Jo 


[ b \ 



+ 

+2 



b2 • n 



(bi • n)(b2 ■ n) 



62 



bi ■ n 



b2 ■ n 



(B21) 



If we choose bi to be parallel to the projection of n into the 
plane of the annulus then bi • n/6 — sin 9vo and b2 ■ n/6 = 0, 
so 



We next average over the velocity distribution of Vm. The 
total rate of energy change is then 



dt 



3^2 J 27r cos" Ovg 



1 + A2 



(1+A")ln(l + A")- A" ^^__2 

A2(l + A2) 

V„,-nln(l + A"^ 



sm tlvo 



Vo A2 



cos 6)14, } /(V)d='V.(B26) 



In general it seems that this equation is not analytically 
solvable, even if /(Vm) is an isotropic Gaussian. However, 
if we are interested in systems where random motions are 
much smaller than the bulk motion (such as galaxy disks), 
then we can approximate /(Vm) = 5(Vm — Vd), where Vd 
is the disk bulk velocity and 5 is the Dirac delta function. 
Note that Vo = Vm — Vm where Vm is the velocity of M. 
For this case 



dt 



,,3,2 I 27rcos Bvo 



1+A2 

(1+A")ln(l + A")- A" 
A2(l + A2) 



(1 — cos 9vo) 



-2nva cos 9^^ ^"'•^^a^ ^ cos 9vo } , (B27) 



cos" 9i,d(l> = TT sin" 9vo ■ 



(B22) 



Using a similiar approach, it is simple to show that 
^ cos9bd<j> = 0. Thus, the change in energy becomes 



ASfl = — {27r AK^i I cos" 9vo + 7rAV;"x sin" 9vo 

+47rVm • nAVmii cos 9va } • 
Substituting eqns. (B15) and (B16) we find 









'-'max 



+47rA' 



bL 



1 +A" 



cos 9vn 



sm 9vo 



(B23) 



where cos 9vo ~ «d cos 9^^ —Vm cos 9^^ . Here «d = Vd/Vb and 
9v^ is the angle between n and Vd, with similar definitions 
for Vm and 9^^^ . The efficiency of energy transfer to direction 
n is then easily found by dividing the above by the same 
expression summed over three orthogonal directions (taking 
one of these to be parallel to Vo simplifies the summation) : 



2 cos" 9v, (1+A")ln(l + A")-A" 2 . . 

H . o^-, — — ( 1 — cos 9va ) 



1 + A2 



A2(l + A2) 



+2 



-2vd cos 9y^ ^^^^j^2^ ^ 
(1 + A")ln(l + A")- A" 



}/{r 



+ A2 



-2vd cos 9v 



A2(l + A2) 
ln(l + A") 



A2 



(B28) 



Vo 



1 +A' 



62 



cos UVr 
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To find the total energy change we multiply by the flux of 
particles passing through the annulus, nVobdb, and integrate 
over 6 from to femax- This gives 



di 



1 + A2 



(1 + A")ln(l + A")-A" . 2, 
+7r^ .o^, — r^TT^ sm 9vo 



'2-K 



A2(l + A2) 
Vm • n ln(l + A") 



Vo 



A2 



cos 9 



Vo 



(B25) 



We are interested specifically in the vertical velocity 
dispersion of a galactic disk. In this case Vd lies in the disk 
plane, while n is perpendicular to that plane. Consequently 
cos 9v^ = and the above expression simplifies to 



2 cos" 9vo _^ (1 + A") ln(l + A 



1 + A2 
2 

1 + A2 



A2(l + A2) 
(1 + A")ln(l + A") 



2. 

(1 — cos 9va 



A2(l + A2) 



(B29) 



This expression is then used in eqn. (2) to calculate the 
vertical heating rate of galaxy disks. Note that < < 1, 
as expected for an efHciency factor. 
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APPENDIX C: DISK SURFACE ENERGY 
DENSITIES 

We here derive expressions for the different contributions to 
the surface energy density of the disk. These are used in 
§2.2.2. 

We assume a disk with a density structure 



p4R,z) = E(i?) 



sech^(z/ij"d) 

Mi 



(Cl) 



with Hd constant with radius. Assuming that the disk is 
thin, Hd <^ R, the potential of the disk can be found by 
approximating the density distribution as a set of infinite, 
homogeneous planes, such that 

/oc 
2TvGpdiR,z)\z - z'\dz' + MR,0) 
-oo 

= 2irGY;{R)H4\ncosh{z/Hd) + ln2] + MH(^) 

Close to the disk plane, the z-component of the force due to 
the spherical halo plus bulge is 



_ GA'UR) 
- ^3 z, 

hence the potential due to these components is 

4'h{R,z) 



(C3) 



GAUR) / , / ^ , 

-z dz + 4>h(R, 0) 



i?3 



2^3 — z + <?!>h(i?,0). 



(C4) 



Referencing all energies to 2: = 0, we can neglect the final 
terms in the above equations. We can now calculate the dif- 
ferent contributions to the disk vertical energy per unit area, 
where necessary using the thin disk approximation Hd <C R- 
The gravitational self-energy of the disk is then 

1 /■°° 

■wdd{R) = - / (l)d{R,z)pd{R,z)dz 

■J —00 

= TvGEiRfHd. (C5) 
The disk-halo gravitational potential energy is 



(C6) 



wdb{R) = / cf>h{R,z)pdiR,z)dz 

■J —00 

and the kinetic energy of the disk is 
U{R) = \y.{R)gI{R). 



(C7) 



numbers of particles. In both cases, the position and veloc- 
ity of the satellite are well converged up until the very final 
stages of the satellite's life (at which point it becomes diffi- 
cult to measure these quantities accurately from the simu- 
lation anyway). The higher resolution simulations lose mass 
from the satellite somewhat more rapidly at late times, but 
the differences are minor and the mass loss rate is well de- 
termined by the simulations. 

The convergence behaviour seems poorer for the change 
in vertical kinetic energy ATz (in which the energy of the 
unperturbed disk from model GISO has been subtracted 
off). For the model G1S5, the difference in the final ATz 
of (0.02 - 0.03) X lO^^Mekm^s"^ between the highest and 
lowest resolution runs could result mostly from the error in 
the energy of the unperturbed disk that is subtracted off, 
since the variation in this value between different realiza- 
tions is around (0.02 - 0.03) x lO^'^Mokm^s"^ in the low 
resolution case. However, for model GlSl the differences in 
ATz between the high and low resolution runs are much big- 
ger than can be explained by errors in the subtraction of the 
unperturbed disk contribution. In this case, the behaviour 
of ATz is not even monotonic as the number of particles 
is increased. We have investigated this further by repeating 
some of the lowest resolution simulations using a different 
sequence of random numbers in generating the initial con- 
ditions. We find that this leads to significant variations in 
ATz, comparable to those seen between the lowest resolu- 
tion and higher resolution simulations. It therefore seems 
that the amount of disk heating by satellites is sensitive 
to stochastic variations in the initial conditions, over and 
above the two-body relaxation which heats the unperturbed 
disk. Comparing the highest and lowest resolution runs, we 
find that the error in the low-resolution estimate of ATz is 
~ 30% for model G1S5, but ~ 100% for model GlSl. The 
convergence of ATz with increasing particle number thus 
seems to depend on the orbital properties of the satellite, 
with different numbers of particles being required to achieve 
the same degree of convergence in different cases. A more 
comprehensive study of convergence in a variety of models 
will be required to address this question fully. 



APPENDIX E: ORDER OF MAGNITUDE 
ESTIMATES OF DISK THICKNESSES 

In this appendix we make order of magnitude estimates for 
the thicknesses of disks resulting from heating by satellites 
and by stars scattering from giant molecular clouds. 



APPENDIX D: CONVERGENCE TESTS ON 
N-BODY SIMULATIONS 

We repeated all the simulations of models GlSl to G1S15 
with one half and one quarter the number of particles, la- 
belling these runs G1S1^''^ GlSl^'"' etc. (In each case, we 
scaled the softening in proportion to the mean interparticle 
separation.) Figures Dl and D2 compare the results for two 
representative models (GlSl and G1S5) run with different 



El Heating by Satellites 

We assume a distribution of satellite halo masses dN jdM = 
AM~^'^, where A is a constant, with a maximum mass of 
/max = 0.01 of the total halo mass and making up a frac- 
tion /mass = 0.1 of the total halo mass (consistent with 
numerical simulations; Springel et al. 2001). Then, 

A = 0.3/ma..A/^'l//L'x. (El) 
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Figure Dl. Convergence tests for model GlSl. Results are shown for the standard simulation (circles), model GlSl^/'^ (which has half 
the number of particles of GlSl; crosses) and model GlSl^^^ (which has one quarter the number of particles of GlSl; triangles). Points 
are connected by dotted lines to guide the eye only — lines are not intended as a realistic interpolation of the points. Top left-hand panel: 
The orbital position of the satellite as a function of time. Top right-hand panel: The orbital velocity of the satellite as a function of 
time. Lower left-hand panel: The remaining bound mass of the satellite as a function of time. Lower right-hand panel: The change in 
the vertical component of the disk kinetic energy due to heating by the satellite as a function of time. 



Taking the Chandrasekhar formula for dynamical friction 
(eqn. B6), the power deposited by a satellite is, to order of 
magnitude 



P = 47r InApd- 



(E2) 



Assuming that heating is effective primarily within one scale 
height of the disk, then heating occurs overly a fraction of 
roughly Hd/2nroYh of the satellite's orbit, where rorb is the 
orbital radius. Allowing for an efficiency of vertical heating 
ez and integrating over masses, the total energy transfer in 



time t is given by 

p , Afd Rdh 
^ = 2t— InA „o, ez 



/maxA/halo 



= 2^1nA-^e.:^0.3/nrassMh'alo/max. (E3) 

Equating to the energy of the disk, as given by eqn. (8) 
after ignoring the contribution proportional to as this is 
typically small, and solving for h 



- 37r"-^^'"^==^"""'' ^Ki?drorbMd ' 



(E4) 
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Figure D2. Convergence tests for model G1S5. Results are shown for the standard simulation (circles), model G1S5^/^ (which has half 
the number of particles of G1S5; crosses) and model G1S5'^'''* (which has one quarter the number of particles of G1S5; triangles). The 
different panels are as for Figure Dl. 
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E2 Heating by Scatterings from Giant Molecular 
Clouds 

We begin with eq. (10), which we approximate to order of 
magnitude (replacing Ae^/At with ez/i etc.) as 



E : 



u\nKal{l3)K,{p)t. 



(E6) 



Inserting eqn. (7) to eliminate (we again ignore the con- 
tribution proportional to h?) this reduces to 

2 GA//dA'/c Sc 



E 



■\nKvat{l3)K,{P)t. 



(E7) 



TT Rih Ed 

Equating to the energy of the disk and solving for h results 



on y\.r 

where / = Ec/Ed. This can be expressed as 
h = 7.2 X 10" 



(E8) 
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